program program_one_turn_map
use madx_ptc_module
use pointer_lattice
use c_TPSA
implicit none


!@ <br> INTERFACE </br>
!@ <br> SUBROUTINE BUILD_LATTICE_ALS(ALS,MIS) </br>
  interface
       subroutine build_lattice_als(ALS,MIS)
         use madx_ptc_module
         use pointer_lattice
         implicit none
         type(layout), target :: ALS
         logical(lp) mis
       end subroutine build_lattice_als
  end interface
!@ <br> END SUBROUTINE BUILD_LATTICE_ALS(ALS,MIS) </br>
!@ <br> END INTERFACE </br>


type(layout), pointer:: ALS
real(dp) prec,closed_orbit(6),mat(6,6)
type(internal_state),target :: state
logical(lp) :: mis=.true. 
integer i,map_order,mf

!!!!!!!!!!!!!!!!!!!!
type(real_8) y(6)
type(damap) mr
type(dragtfinn) dfr
!!!!!!!!!!!!!!!
type(c_taylor) ct,dct,ict
type(c_normal_form) normal_form
type(c_vector_field) df
type(c_factored_lie) fn
type(c_damap)  one_turn_map, Id,m1
integer pow(lnv)
prec=1.d-6 ! for printing
longprint=.false. 
 
call ptc_ini_no_append
call append_empty_layout(m_u)
ALS=>m_u%start

call build_lattice_als(ALS,mis) 

state=delta0

map_order=3
call init_all(state,map_order,0)
write(6,*) " "
call alloc(one_turn_map,id,m1)
call alloc(y) 
call alloc(normal_form) 
call alloc(ct,dct,ict) 
call alloc(df)
call alloc(mr)
call alloc(dfr)
call alloc(fn)

closed_orbit=0.d0
call find_orbit_x(als,closed_orbit(1:6),STATE,1.e-5_dp,fibre1=1)     ! (1)

id=1   ! map is set to identity                                      ! (2)

! map is added to closed orbit and put into the 6 polymorphs
y(1:6)=closed_orbit(1:6)+id                                          ! (3)

call propagate(als,y(1:6),state,fibre1=1)                            ! (4)

one_turn_map=y(1:6) ! Six polymorphs are promoted to Taylor maps     ! (5)
mr=one_turn_map
closed_orbit=y                                                       ! (6)




call  c_normal(one_turn_map,normal_form)                             ! (8)


!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!@3  <a name="CMONO.2"></a>
!!! ct= (3+5i)*z_2
write(6,*) "ct= (3+5i)*z_2"
ct=(3.d0+5*i_).cmono.2
call print(ct,6)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!

write(6,*) "ct= (2+6i)*z_2*z_3^2"
POW=0; POW(2)=1;POW(3)=2;
ct=(2.d0+6*i_).cmono.POW
call print(ct,6)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!

write(6,*) "ct= (2+6i)*z_2*z_3^2"
ct=(2.d0+6*i_).cmono.'012'
call print(ct,6)


!!! ct= (3.d0)*z_2
write(6,*) "ct= (3)*z_2"
ct=(3.d0).cmono.2
call print(ct,6)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!

write(6,*) "ct= (2.d0)*z_2*z_3^2"
POW=0; POW(2)=1;POW(3)=2;
ct=(2.d0).cmono.POW
call print(ct,6)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!

write(6,*) "ct= (2.d0)*z_2*z_3^2"
ct=(2.d0).cmono.'012'
call print(ct,6)
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
write(6,*) " "
write(6,*) "ct= (2.d0)*z_2*z_3^2 + (2.d0+i)*z_2*^2*z_3^1"
ct=((2.d0).cmono.'012')+((2.d0+i_).cmono.'021')
call print(ct,6)
dct=ct.d.3
write(6,*) "dct = ct.d.3 = d ct / dz_3"
call print(dct,6)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
write(6,*) " "
ict=dct.i.3
write(6,*) "ict = \int dct dz_3 "
call print(ict,6)


!!!!!!!!!!!!!!!!!!  factored maps with nonlinear part as one Lie exponent !!!!!!!!!!!!!!!!!!!!!!!
!type c_vector_field
call kanalnummer(mf,"One_Lie_factored.txt")
call c_factor_map(one_turn_map,m1,df,1) 
write(mf,*) " M= df o m1 "
call print(m1,mf,prec)
write(mf,*) "   df = exp(f.grad) Id"
write(mf,*) "    f is printed "
call print(df,mf,prec)
dfr=mr
call print(dfr,mf)
call c_factor_map(one_turn_map,m1,df,-1) 
write(mf,*) "  M= m1 o df"
call print(m1,mf,prec)
write(mf,*) "   df = exp(f.grad) Id"
write(mf,*) "    f is printed "
call print(df,mf,prec)
close(mf)

!!!  factored maps as produced by normal forms:    multiples products of nonlinear vector fields !!!
!!!! c_factored_lie
!!! this is done naturally by the normal form
call kanalnummer(mf,"Dragt_Finn_factored.txt")
write(mf,*) " the nonlinear part canonical transformation A of the normal form "
write(mf,*) " factorised as : exp(normal_form%g%f(no).grad) ... exp(normal_form%g%f(1).grad) "
fn=normal_form%g
call print(fn,mf,prec)
close(mf)

!!!!!!!!!!!!!!!! exponential operators !!!!!!!!!!!!!!!
!!!!  M= m1 o df"   implies     M = exp(f . grad) m1
!!!!   let us try that  
call kanalnummer(mf,"exponential_operators.txt")

write(mf,*) " Testing the one Lie factorization "
id=texp(df,m1)  !!! notice the reverse order of Lie operators
id=id**(-1)*one_turn_map

call print(id,mf,prec)
!!!!!!!!!!!!!!!!   
!!!            n%a_t=n%a1*n%a2*from_phasor()*texp(n%g)*from_phasor(-1)
!!!  texp(n%g) is same as  texp(n%g,Identity)

!!! Let us try using Lie operators
write(mf,*) " Testing the factored Lie factorization "
id=normal_form%a1*normal_form%a2*from_phasor()
id=texp(normal_form%g,id)  !!! notice the reverse order of Lie operators
id=(id*from_phasor(-1))**(-1)

id=id*normal_form%a_t
call print(id,mf,prec)
close(mf)


call ptc_end

end program program_one_turn_map
